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We address qqQQ exotic tetraquark bound states and resonances with a fully unitarized and 
microscopic quark model. We propose a triple string flip-flop potential, inspired in lattice QCD 
tetraquarks static potentials and fluxtubes, combining meson-meson and tetraquark potentials. Our 
potential goes up to the color excited potential, but neglects spin-tensor potentials. To search for 
bound states and resonances, we first solve the two-body mesonic problem. Then we develop fully 
unitary techniques to address the four-body tetraquark problem. We fold the four-body Shcrodinger 
equation with the mesonic wavefunctions, transforming it into a two-body meson-meson problem 
with non-local potentials. We find bound states for some quark masses numbers, including the one 
reported in lattice QCD. Moreover, we also find resonances and calculate their masses and widths, 
by computing the T matrix and finding it’s pole positions in the complex energy plane, for some 
quantum numbers. 

However a detailed analysis of the quantum numbers where binding exists shows a discrepancy 
with recent lattice QCD results for the Ubb tetraquark bound states. We conclude that the string 
flip-flop models need further improvement. 


I. INTRODUCTION 

A long standing problem of QCD is the search for lo¬ 
calized exotic states [1] and the corresponding decay to 
the hadron-hadron continuum. There is no QCD theorem 
preventing the existence of exotic hadrons, say two-gluon 
glueballs, hybrids, tetraquarks, pentaquarks, three-gluon 
glueballs, hexaquarks, etc, and the scientific community 
continues to search for clear exotic candidates. How¬ 
ever, this problem turned out to be much harder than 
expected. 

In this work, we develop fully unitary techniques to 
solve some of the theoretical problems of multiquarks, 
and predict multiquark bound states and resonances. We 
continue a previous unitary study of tetraquarks with a 
simplified two-variable toy model [2], now fully solving 
the tetraquark problem, of the Schrodinger equation for 
two quarks and two antiquarks. We specialize in systems 
which are clearly exotic tetraquarks, who cannot have a 
significant mesonic quark-antiquark component, i e where 
the quantum numbers, or the S-matrix pole and decay 
amplitudes clearly show it is a tetraquark. 

In particular, as a benchmark, we study in detail the 
light-light-antiheavy-antiheavy systems who are expected 
to produce tetraquarks. From basic principles of QCD, 
it is clear a system with two light quarks and two heavy 
antiquarks, for instance with flavour udbb, should form 
a boundstate if the antiquarks are heavy enough [3-13]. 
To understand why binding should occur it is convenient 
to use the Born-Oppenheimer [14] perspective, where the 
wavefunction of the two heavy antiquarks is determined 
considering an effective potential integrating the light 
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Figure 1: Screening in the light-light 
-antistatic-antistatic tetraquark [15]. 


quark coordinates. At very short bb separations r, the b 
quarks interact with a perturbative one-gluon-exchange 
Coulomb potential, while at large separations the light 
quarks totally screen the interaction and the four quarks 
form two rather weakly interacting B mesons as illus¬ 
trated in Fig. 1 . Thus a screened Coulomb potential is 
expected. This potential clearly produces a boundstate, 
providing the antiquarks bb are heavy enough. 

We leave the brief review of the tetraquarks in the 
literature for Section II, including the searches for 
tetraquarks in experiments, in dynamical lattice QCD 
and in quark models. We also discuss searches combin¬ 
ing semi-dynamical lattice QCD with quantum mechani¬ 
cal techniques. In Section III we describe our triple string 
flip-flop potential, where our system is open to the contin¬ 
uum in the inter meson-meson directions, and is confined 
in the tetraquark direction. Our string flip-flop potential 
also includes the first color excitation. In Section IV, we 
address the meson-meson scattering problem, occurring 
when we solve the Schrodinger equation. We detail our 
numerical techniques, utilized to solve both our bound 
state and scattering problems. In Section V we show 
our results with the state of the art triple string flip-flop 
potential, exhibiting tetraquark bound states, and res¬ 
onances. For the resonances, with our fully unitarized 
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computation, we calculate the pole position and thus 
finding their decay widths. We also consider, to com¬ 
pare with our full computation, simplified potentials. In 
Section VI we compare our results to the existing lattice 
QCD results for the light-light-antiheavy-antiheavy sys¬ 
tem and conclude. We find excessive binding, concluding 
the flip-fop potentials need further improvement. Never¬ 
theless our technique can be applied to other potentials 
who may be developed in the future. 


II. BRIEF REVIEW OF TETRAQUARKS 
A. The experimental search for tetraquarks 
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This is a very difficult problem experimentally, since 
exotic candidates are resonances immersed in the excited 
hadron spectra, and moreover they usually decay to sev¬ 
eral hadrons. 

Recently an experimental article [16] was published 
indicating that the 7ri(1600) observation of a resonance 
with mass M = 1660(74) MeV , and width r= 269 (85) 
MeV, in diffractive dissociation of 190 GeV/c tt— into 
7r“7r“7r“'", has exotic = 1 '' quantum numbers, 

consistent with an hybrid meson or a tetraquark. 

Moreover, the existence of tetraquarks has been ad¬ 
vanced by the experimental collaborations at the charm 
and bottom factories to interpret the new X, Y, Z hid¬ 
den charm or bottom mesons [17]. Their mass and decay 
products mark them as charmonia-like resonances but 
their masses do not fit into the quark model spectrum of 
quark-antiquark mesons [18]. This class of tetraquarks is 
related to the double-heavy tetraquarks we address here. 

In particular, the charged and Z^ are crypto¬ 
exotic, but technically they can be regarded as essen¬ 
tially exotic tetraquarks if we neglect cc or bb annihila¬ 
tion. There are two Z^ observed only by the collabora¬ 
tion BELLE at KEK [19], slightly above B B* and B* B* 
thresholds, the Y;i(10610)“'' and Zb(10650)+. Their na¬ 
ture is possibly different from the two Yc(3940)^and 
Yc(4430)^, whose mass is well above DD threshold [20]. 
The Z^ has received a series of experimental observa¬ 
tions by the BELLE collaboration [21, 22], the Cleo-C 
collaboration [23], the BESIII collaboration [24-28] and 
the LHCb collaboration [29]. This family is possibly re¬ 
lated to the closed-charm pentaquark recently observed 
at LHCb [30]. 

However, to establish a new resonance it is necessary 
to study with an accurate level of confidence all its prop¬ 
erties, including its mass and width as determined by its 
5'-matrix pole and all relevant partial decay widths. Pos¬ 
sibly we need more data and more extensive data analysis 
to be able to absolutely confirm exotics [31, 32]. 

Notice, using very approximate Resonant Group 
Method calculations, in 2008, we predicted [33] a partial 
decay width to tt J/ip of the Yc(4430)“ consistent with 
the experimental value. But we opt to leave the detailed 
study of this second class of double-heavy tetraquarks. 


Figure 2: Surface plot of the static tetraquark flux tube 
computed in lattice QCD [34, 35] . 


with closed charm or closed bottom tetraquarks, for fu¬ 
ture studies. We first want to address more constrained 
systems with identical fermions, where the Pauli symme¬ 
try applies. 


B. The lattice QCD search for tetraquarks 

In lattice QCD, the study of exotics is presently even 
harder than in the laboratory, since the techniques and 
computer facilities necessary to study of resonances with 
many decay channels remain under development. 

Thus lattice QCD started by searching for the expected 
boundstate in light-light-antiheavy-antiheavy channels 
[36, 37]. Using dynamical quarks, the only heavy quark 
presently accessible to Lattice QCD simulations is the 
charm quark. No evidence for boundstates in this family 
of tetraquarks, say for a udbb was found. 

Lattice QCD also searched for evidence of a large 
tetraquark component in closed bottom the Zc(3940)“ 
candidate [38, 39]. The difficulty of the study of the Z ~, 
a resonance well above threshold, is due to its many two- 
meson coupled channels. The authors considered 22 two- 
meson channels, corresponding to lattice QCD interpola¬ 
tors . In addition they considered 4 tetraquarks 

channels, corresponding to diquark-antidiquark interpo¬ 
lators with flavour and color [cm ]3 [cdjg. An evidence for 
the tetraquark resonance candidate was investigated in 
the full coupled correlator matrix of hadron operators. 

Finally, after switching on and off the 4 tetraquark- 
like channels, the authors [38, 39] found no significant 
deviation in the 13 lowest channels, who span the energy 
range from the lowest threshold to the Yc(3940)“ can¬ 
didate. Thus, the authors concluded there is no robust 
evidence of a Z^ tetraquark resonance. 

However, the direct proof for, or against, a tetraquark 
resonance in lattice QCD would require the study of the 
S-matrix. The technique to perform phase shift analysis 
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in lattice QCD exists [40, 41]. From the phase shift anal¬ 
ysis, inasmuch as with experimental data, the poles of 
the S-matrix can be extracted. But phase shift analysis 
of the tetraquark has not yet been done with lattice 
QCD data. For an absolute evidence, the different par¬ 
tial decay widths should be computed as well in lattice 
QCD. 

With present computers only resonances with just ^ 
1 open decay channel have been studied in lattice QCD, 
with sufficient detail [42]. The method of extracting the 
phase shifts from the spectrum of harmonic waves in a 
box has been extended to inelastic (more than one open 
channel) coupled channels [43]. However, tetraquarks 
are excited resonances, decaying into many channels, ^ 
30 for the last experimental Z~ candidates. This is 
presently unattainable by present computers and codes, 
but we expect Lattice QCD to be on the way to, in the 
future, reach the level of experimental data analysis. 



C. Quark model approaches to exotic tetraquarks 

A detailed theoretical understanding of the properties 
of exotic hadrons is important to support the experimen¬ 
tal and lattice QCD searches of exotics. Several theoret¬ 
ical problems need to be solved before addressing exotics 
hadrons. 


1. Early quark models 

Already at the onset of QCD, the bag model predicted 
many tetraquarks [1[. However, as soon as lattice QCD 
was able to compute static quark potentials and color 
electromagnetic fields , it was realized that quark con¬ 
finement was not bag-like, but string-like, due to color 
flux tubes. 

Inspired in lattice QCD linear confining potentials, the 
relativized quark model potential was developed [18, 44], 
after the authors fitted the spectra of all known hadrons 
in the 80’s. Notice that a correctly calibrated quark 
model needs many terms and many parameters, say of 
the order of ~ 20 parameters. 

Nevertheless the relativized quark model still lacks two 
crucial effects, leading with up to 400 MeV deviations, 
chiral symmetry breaking, which has already been for 
instance in the Dyson-Schwinger approach [45-47], and 
coupled channels / unquenching which has already been 
include for instance in effective meson or baryon models 
[48[. 


2. Extra binding with four-body flux tube potentials 

Moreover, since tetraquarks are always open to decays 
into a meson-meson pair, tetraquark resonances or bound 
states may only exist if a mechanism exists to provide 


Figure 3: Triple string flip-flop Potential potential. 
While the previous string flip-flop potentials choose the 
minimum of two different meson pair potentials, we 
consider as well the tetraquark potential [2[. 


binding specific to tetraquarks. Here we explore a mech¬ 
anism observed in lattice QCD static potentials: the con¬ 
fining four body potential [49, 50], produced by double 
Y or butterfly shaped flux tubes or strings [34, 35, 51], 
depicted in Fig. 2. This mechanism is related to the 
Jaffe-Wilczek model [52-54] who proposed the tetraquark 
would form a diquark-antidiquark system. 

We acknowledge other mechanisms may exist to sup¬ 
port binding, however they are complex to implement 
and it is not clear, from lattice QCD, how they work 
quantitatively. For instance attraction may also be due 
to quark-antiquark annihilation, however it turned out 
to be insufficient to bind a proposed pentaquark [55, 56]. 
Another mechanism is the hyperfine spin-dependent po¬ 
tential utilized in the original bag model [1[, however the 
spin-tensor potentials have only been computed in lattice 
QCD for mesons. For baryons they are model dependent, 
while for in multiquarks the details of the spin-tensor 
potentials remain speculative. Moreover both quark- 
antiquark annihilation and the hyperfine potential are 
also important for chiral symmetry breaking. To avoid 
the complexity of quark-antiquark annihilation, spin ten¬ 
sor quark-quark interactions and chiral symmetry break¬ 
ing, we specialize in a family of tetraquarks where they 
can be neglected. 

Here we consider purely exotic tetraquarks only. In 
contradistinction to crypto-exotics, in pure exotics quark- 
antiquark annihilation does not occur directly. Crypto¬ 
exotic tetraquarks are also less clear in the sense they 
always have a mesonic component, they are never a pure 
exotic. Thus we consider only tetraquarks with exotic 
flavour. Moreover, as a first study, we neglect spin-tensor 
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interactions and chiral symmetry breaking effects. 


3. Solving the Van der Waals force problem with string 
flip-flop potentials 


Clearly, tetraquarks are always coupled to meson- 
meson systems, and we must be able to address correctly 
the meson-meson interactions. 

Notice confining two-body potentials with the SU(3) 
color Casimir invariant Xi ■ Xj suggested by the One- 
Gluon-Exchange potential, and possibly compatible with 
lattcie QCD, lead to a Van der Waals potential, 


Vv: 


an der Waals — 


W(r) 


X T. 


( 1 ) 


where T is a polarization tensor. This would lead to 
an extremely large Van der Waals [57-62] force between 
mesons, or baryons which clearly is not observed experi¬ 
mentally. Thus two-body confinement dominance is ruled 
out for multiquark systems. 

The string fiip-fiop potential for the meson-meson in¬ 
teraction was developed [52, 63-66], to solve the problem 
of the Van der Waals forces produced by the two-body 
confining potentials. 

Traditionally, the string fiip-fiop potential considers 
that the potential is the one that minimizes the energy 
of the possible two different meson-meson configurations, 
say Mi 3 M 24 or M 14 M 23 . This removes the inter-meson 
potential, and thus solves the problem of the Van der 
Waals force. 

Here we upgrade the string fiip-fiop potential, con¬ 
sidering a third possible configuration, the tetraquark 
one, say 7 i 2 , 34 , where the four constituents are linked 
by a connected string [13, 67]. We study whether the 
tetraquark attractive flux tube may induce further bind¬ 
ing of tetraquarks. 

The three configurations differ in the strings linking 
the quarks and antiquarks, this is illustrated in Fig. 
3. When the diquarks qq and qq have small distances, 
the tetraquark configuration minimizes the string energy. 
When the quark-antiquark pairs qq and qq have small 
distances, the meson-meson configuration minimizes the 
string energy. 
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Figure 4: Color of the diquarks determined by the 
quantum numbers of BB tetraquark attractive channels 

[71[. 


With a triple string fiip-fiop potential, bound states 
below the threshold for hadronic coupled channels have 
been found [13, 67-69]. 

On the other hand, the string flip-flip potentials allow 
fully unitarized studies of resonances [65, 66 , 70]. Uti¬ 
lizing analytical calculations with a double fiip-fiop har¬ 
monic oscillator potential [70], and using the resonating 
group method again with a double fiip-fiop confining har¬ 
monic oscillator potential [65, 66 ], resonances and bound 
states have already been predicted. 

Moreover, using the perturbative approximation of the 
resonating group method, a preliminary estimation of the 
partial decay width of the Z(4430)“ resonance was sim¬ 
ilar to the one measured by LHCb [33]. 

These studies suggest tetraquark bound states or reso¬ 
nances are plausible. Fully unitarized techniques adapted 
to state of the art potentials, remain to be applied 
in order to achieve a quantitative theoretical study of 
tetraquark resonances and bound states. 


D. Recent lattice QCD results for the potentials in 
the light-light-antistatic-antistatic system. 


f. Previous quark model results with string flip-flop 
potentials 

Tetraquark resonances are quite subtle, since the Fock 
space of tetraquarks is the same of its decay channels, the 
meson-meson channels, and moreover a potential barrier 
is absent. A priori it is not intuitive whether this system 
may produce resonances or not. 

Nevertheless, there is an argument [52] suggesting mul¬ 
tiquarks with angular excitations may gain a centrifugal 
barrier, leading to narrower decay widths. 


The difficulties of both lattice QCD and quark models 
can be both relaxed when one uses a hybrid approach 
with both numerical and model techniques. 

Recently, the light-antistatic light-antistatic potentials 
have been computed in in lattice QCD [72, 73], as shown 
in Fig. 5 . A static antiquark constitutes a good 
approximation to a spin-averaged b bottom antiquark. 
The light-antistatic light-antistatic potential can then be 
used, with the Born-Oppenheimer approximation [14], as 
the & B — B potential, where we the higher order 1/mb 
terms including the spin-tensor terms are neglected. Ac¬ 
cording to the quantum numbers of the two dynamical 
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(a) scalar isosinglet:« = 0.29 ± 0.03, p = 2.7 ± 1.2. d/a = 4.5 ± 0.5 



(b) vector isotrlplefio = 0.20 ±0.08, d/a = 2.5± 0.7 (p = 2.0 fixed) 



probability to find the b antiquark pair at separation r 



Figure 6 : Probability density of the bb pair determined 
with the Born Oppenheimer approximation [15]. 


light quarks, either attraction or repulsion is found. The 
attraction/repulsion can be understood with the screen¬ 
ing mechanism illustrated in Fig. 1 , and the antistatic- 
antistatic potential is well fitted by the screened Coulomb 
ansatz [15], 



Utilizing the potential of the channel with larger at¬ 
traction, occurring in the lsospin=0 and Spin=0 quark- 
quark system, together with the Born-Oppenheimer ap¬ 


proximation [14] to derive the Schrodinger equation, 
the possible boundstates of the heavy antiquarks are 
then investigated with quantum mechanics techniques. 
Recently, this approach indeed found evidence for a 
tetraquark udbb boundstate [15, 74], while no bound- 
states have been found for states where the heavy quarks 
are cb or cc (consistent with full lattice QCD computa¬ 
tions [36, 37[) or where the the light quarks are ss or 
cc [71]. The bb probability density in the only binding 
channel is shown in Fig. 6 . 

Moreover, lattice QCD finds attraction or repulsion, 
consistently with the screening model. 

• In some channels we find attraction, in others we 
find repulsion. 

• Phenomenologically, due to the short range one- 
gluon exchange proportional to the Gell-Mann A® • 

Casimir operator, the b and b are expected to 
bind only if they are in a triplet 3, not a sextet 6 , 

• thus the light quarks qq should be in a color an¬ 
titriplet 3, as illustrated in Fig. 4, 

• due to the Pauli principle and to the quark model 
phenomenology this is indeed consistent with an 
s-wave, 1=0 and S=0 diquark. 

These results for repulsion, attraction or binding are very 
important, the quark models should comply with them. 


III. POTENTIAL 

A. Groundstate string flip-flop potential 

We know from lattice computations for static quarks 
[34, 35, 49, 50] that the ground state potential for a sys¬ 
tem composed of two quarks and two antiquarks, is well 
fitted by a string flip-flop potential , 

Vpp = min(VMM, Umm', Ur) , (3) 

were, for simplicity, we neglect possible mixings at the 
transition regions. In Eq. (3), Vmm and Vmm' are two 
possible potentials of a pair of mesons, given by the sum 
of the intra-meson potentials UMM(xi) = Um(|xi — xaj)-|- 
Um(|x 2 -X 4 |) and UMM'(xi) = Um(|xi-X 4 I)-t Um(|x 2 - 
X 3 I). The intra-meson potential is well described by the 
funnel potential, 

UM(Ui) = K - — + (TTij , (4) 

Tij 

where Vij = jx^ — Xj| and which includes a constant term 
K, the short range Coulomb potential — - and the long 
range confining potential err. Here we use the indices 1 
and 2, to refer to the quarks while 3 and 4, refer to the an¬ 
tiquarks. Xi are the positions of the quarks/antiquarks. 
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(a) Meson-meson groundstate flux tube 


(b) Tetraquark groundstate flux tube 


Figure 7: Lagrangian density plots for static quark-quark-antiquark-antiquark flux groundstate tubes computed in 
lattice QCD [35]. The density and the coordinates are in lattice spacing units [35]. Three different flux tubes occur 

for different geometries of the four-body system. 
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(a) Meson-meson excited flux tube 


(b) Tetraquark excited flux tube 


Figure 8: Lagrangian density plots for static quark-quark-antiquark-antiquark first excitation flux tubes computed 
in lattice QCD [35]. The density and the coordinates are in lattice spacing units [35]. We show the first excitation of 
the different groundstate flux tubes shown in Fig. 7. In 8a the excitation of a meson-meson flux tube is similar to the 
other meson-meson flux tube. In 8b the excitation of a tetraquark flux tube is similar to a meson-meson flux tube. 


Moreover the tetraquark function Vr is, 

4 

Vr(x,) =2K- Cb— + , (5) 

l=^<J 

where Cij = 1/2 for quark-quark and antiquark- 

antiquark and Cij = 1/4 for quark-antiquark interac¬ 
tions. Lmin is the minimal distance linking the four par¬ 


ticles, 

LminiXi) =ri5+r25 + ^5 6 + J'S 6 + ?'4 6 , (6) 

where 5 and 6 are the indices of the two Fermat-Torricelli- 
Steiner points [13, 67, 75] of the tetraquark. We compute 
the position of these two points with the numerical tech¬ 
nique of Ref. [75]. This potential generalizes the earlier 
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string flip-flop potential models by introducing a third 
possible branch in the potential where the four particles 
are all linked by a confining string. 

Our main approximation consists in neglecting the spin 
tensor potentials, say the hyperflne, spin-orbit and ten¬ 
sor potentials. These potentials have only been studied 
in detail in the two-body mesonic quark-antiquark sys¬ 
tem, including in lattice QCD computations [76]. In the 
three quark baryonic system, the spin tensor potentials 
have been applied in models, however they remain to 
be computed in lattice QCD. In four-body tetraquarks 
systems, we still ignore the spin-tensor potentials. More¬ 
over, the study of tetraquarks with a four-body confining 
potential remains an open problem, and our present aim 
is to solve it before addressing quark spin interactions. 
Thus our tetraquark masses or poles should be seen as 
spin-averaged results. 


It is useful to introduce the contravariant basis to \Ci) 
and |C//). The color space metric matrix in this basis is. 


9 = 


(I 1/3 


and it’s inverse 


9 = 


1/3 

1; 

9/8 

- 3/8 

-3/8 

9/8 

while 

5-1 is 


( 10 ) 


( 11 ) 


ing the metric matrices, we define the contravariant basis 
\C^) = g^^\CB), with the property {C^\Cb) = Sq- These 
contravariant basis is explicitly defined by 


\C^) = g|C/) - 8’ 

lO = -^|C/) + ^|C//) ■ (12) 


B. Color states 

However, this potential is not yet sufficiently defined 
to completely address tetraquarks. It has been pointed 
out [77-81] that we also need a potential for color excited 
states, otherwise there would be solutions which wouldn’t 
be color singlets. 

With two quarks and two antiquarks, two different 
linearly independent color singlets can be constructed. 
Physically, rather than choosing a basis composed by 
singlet-singlet and octet-octet, or including the triplet- 
antitriplet, it is more convenient to have a color basis cor¬ 
responding to the two possible asymptotic meson-meson 
systems MM and MM'. In the color system denoted 
with I the pairs Q 1 Q 3 and Q2Q4 form color singlets, and 
in the color system denoted with II, the pairs Q 1 Q 4 and 
Q2Q3 are color singlets. These two color states are given 
by: 

1^/) ^^CiC3^C2C4 IC1C2C3C4) 

1^//) — ^^CiC4^C2C3 IC 1 C 2 C 3 C 4 ) (7) 

However, these two states are not orthogonal, (C/|C//) = 
|. For instance an orthogonal basis could be constructed 
by considering the antisymmetric A and symmetric S 
color combinations of these two states, 

l“4) = ~ 

\S) = ^{\Ci) + \Cu)) ( 8 ) 

and inversely, 

|C//) = yi|5)-^|M> (9) 


C. Completing the potential 

Since we have two independent color singlets, the po¬ 
tential must be given by a 2 x 2 matrix. The lowest 
eigenvalue of the matrix should correspond to VpF- 

The corresponding eigenvectors |mo) are the following, 

• when VpF = Vmm, the corresponding eigenvector 
is |mo) = \Ci) , 

• Vff = Vmm' the corresponding eigenvector is 
\uo) = \Cn) , 

• when Vff = Vt, we have jwo) = |M). 

The eigenvectors for the highest eigenvalue jrti) must 
be orthogonal to [mq) so that the potential matrix is her- 
mitian. So, in the three cases, we have the following 
eigenvectors, 

• when = Vmm, |mi) = \Ci) = 

• when Vff = Vmm', \ui) = \Cn) = 

• when Vff = Vt, |mi) = |5) 

To construct the potential matrix, we only need to 
know the highest eigenvalue. In Refs. [35, 80], in the 
transition region, there seems to be an exchange between 
the ground and the first excited states. Namely, in re¬ 
gion A (where Va is the ground state) close to the tran¬ 
sition to region B (where Vb is the ground state), the 
ground state’s potential is vq = Va and the first excited 
state Vi = Vb- When we enter region B, then we have 
vq = Vb and vi = Va- This way, the first excited state 
should be the second lowest of the three potentials close 
to the transitions. We assume, that this result is valid in 
general, i e 

VpF = min[max(VMM, Vmm')> inax(t4fM, Vt), 

max(VMM',hT)] ■ (13) 



D. Other Models 

For comparison, we also use two other models. One of 
them is similar to this color structure dependent triple 
flip-flop model, but with the tetraquark sector removed 
from the potential. So, the ground state’s potential is 
just 

uo = min(V 7 , Vf/) , (14) 

and the excited state is, 

ui = max(V 7 , V//) , (15) 

with the corresponding eigenvectors, being constructed 
the same way as before. This potential is given in the 
matrix form on the non-orthogonal basis formed by \Ci) 
and \Cii) by 

T/ = + I ^ii) I min(V/, Vu) \ 

imin(yf,Vf/) IVii + ) 

(16) 

A third model we also compare with, is the colorless 
double Flip-flop. In this model, the potential is simply 
given by 

V = min{Vi,Vii) . (17) 

This potential does not depend of the color degrees of 
freedom of the system and so is a non-physical one. 
However, this model and it’s extension to include the 
tetraquark sector have been used by several authors 
[65, 66, 70]. It can be interpreted as the limit of a color 
structure dependent model, when the two eigenstates of 
the potential are degenerate (uq = ^i)- 


IV. OUR UNITARY TECHNIQUE TO SOLVE 
THE MESON-MESON SCATTERING AND FIND 
TETRAQUARK BOUND-STATES AND 
RESONANCES 


A. Meson-meson scattering equation 


Let us start with the microscopic Hamiltonian, 

H = fQ + VQ, (18) 


where the Q subscript means that we are dealing with 
the kinetic and potential energy of the quarks and not 
of the mesons. Folding the operators in the color states 
|C/) and |C//), we get the Schrodinger equation: 

{9abTq + = EgAB-^^ . (19) 


We want to study meson-meson scattering because 
mesons and not quarks are the asymptotic states of the 
theory. For this, we need the meson kinetic energy opera¬ 
tor to show up in the Schrodinger equation. The mesons 
kinetic energy is given by 


Tm 


Tq + Vmm 

0 


0 ") . 

Tq + Vmm' J 


W now need to opt for a decomposition of the potential 
and kinetic energies in the Hamiltonian. 

However, due to the non-orthogonality of our color ba¬ 
sis, it is not trivial to conveniently decompose the Hamil¬ 
tonian. A first guess would be to consider the meson 
potential energy as. 


wM _ 
^AB — 


Vi I Viii 
Vu I Vu u 


f 1 1/3 \ f Vmm 0 

I 1/3 1 [ 0 Vmm' 


Vi I — Vmm Vi n — | Vmm> 
Vii I — \ Vmm Viiii — Vmm' 


( 20 ) 


however this turns out to be a non-Hermitean operator, 
and clearly this is not a good decomposition of the Hamil¬ 
tonian. Another possible decomposition would be to add 
the potentials Vmm and Vmm' not to TqI but to gTq, 


( Tq + Vmm ^Tq \ 
\ ^Tq Tq + Vmm' J 


however we do not want Tq to appear in the equa¬ 
tion. We could as well define Tq = T^ — which 
is convenient since in one case the operator is applied 
to the two components of the wave-function. However, 
numerically it could be difficult to verify the equality 
T/f — VIj = T//i — . Moreover, a simple reduced 

mass correction or the use of a relativistic quark energy 
destroys this equality. 

After considering several decompositions and testing 
them, we Anally solve all these problems by considering 
the decomposition. 


Tmm+Tmm' 
Tc = I - 6 

V Tmm' 


and 

Vs = 


Vi I — Vmm 

T/ Vmm+Vmm' 

yj 11 Q 


Vi II- 


^MM + yMM' 


Vii II — Vmm' 


( 21 ) 


( 22 ) 


Notice, both operators are now Hermitean. We have the 
Schrodinger equation. 


+ Vab'^^ = EgAB^^ ■ 
The components are expanded as, 

= X! (^13, (ri324) , 

i 

i 

where we define the Jacobi coordinates, 

Pi3 = Xi - X3 , 

P24 = X 2 — X 4 , 

TOiXi + 7743X3 7772X2 + 7774X4 

ri324 — -,-,- 

7771 + W 3 7772 + 7774 

Pl4 = Xi — X 4 , 

P23 = X 2 — X 3 , 

777iXi + 7774 X 4 7772 X 2 J- 7773 X 3 

ri423 — - 1 -;- 

777i J- 7774 7772 + 7773 


(23) 


(24) 


(25) 
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The functions are eigenfunctions of Hamiltonian of 
the two non-interacting mesons, 

{Ta + VA)<i>f = . (26) 

They have the property, 

lim $f(pi,p 2 )= lim ^f{pi,p 2 ) = 0 , (27) 

Pi —¥00 p 2—¥00 

because the quarks in each meson are confined. 

Substituting Eq. (24) into Eq. (23) and integrating the 
confined degrees of freedom, (i. e. the pi coordinates), 
we obtain the multichannel equation, 

. (28) 


Here, we employ Greek indices to denote both the color 
structure and internal quantum numbers. 

The operators Vap and ga/j are local between 

states with the same color structure but are non-local 
between states with different color structures Let us now 
analyze in detail the non-local operators. Consider the 
non-diagonal term of the g matrix, folding it with and 



^'*(Pl3,P24)*V'^ (1-1324)* 

‘J’j’^(Pl4, P23)V'^^( 1 - 1423 ) ■ (29) 

Changing the integration variables to ri 324 , ri 423 and 


miXi -I- 7712X2 TO 3 X 3 -I- 7774 X 4 

1-1234 — -,-,- 

7771 + 1172 III3 + W4 


(30) 


in order to isolate functions of ri 324 and ri 423 , the inte¬ 
gral becomes 


J rf^ri324fi^ri423 Ipi (ri324)*7ij (1-1324, l-i423)V'j^(l-1423) ■ 


The function 7ij(ri324, ri423) is given by. 


7y- = A 


<i^ri234 (Pl3,P24)*‘kj^(Pl4,P23) , 


(31) 


where. 


/ Mi2M^4^Mi^M2aMhM2Z \ ^ 

\ 2r77i77l27773r774M2 / 


(32) 


with Mij = rrii + rrij and M = 777 ^. In this way, 

knowing that g acts on a function as, 

5a/3V'/3 = J C^^l-/35a/3(l-a,l-^)l/’'^(l-^) , (33) 

we have, 

ff/i,/j(ri324,1-1324) = ^ij '^^(ri324 ~ ri324) , 
3/i,//j(ri324, 1-1423) = g7ij (1-1324, 1-4423) , 

3//i,/j(ri423, 1-4324) = g7ii(l-1324, 1-1423) , 
ff//*,7/7 (1-1324, r'1423) = '^^(1-1423 “ l-'l423) ■ (34) 


The potential V has a similar structure, 

VA^,AJirA,rA) = V^^irA)S^{rA - y'a) , 
j ( 1 - 1324 , 1 - 4423 ) = Cj ( 1 - 4324 , 1 - 4423 ) > 

L 77 i,/j(ri 324 , 1 - 4324 ) = l'4i(l-1324, 1 - 1423 ) , (35) 

with 

^ J '^^Pl3'^^P24 4>f(pi3,P24)*(V/y - yf)$j(pi3,p24) , 

= y C?^Pl4C?^P23 4 >[(Pi4, P23)*(L//.// — 47/)4>}(P14, P 23 ) , 
Vij = A y rf^ri234 

4^f(Pl3,P24)*(V/,7/- g ^^ )4>j^(Pl4,P23) • (36) 

Note that, because of Eq. (27), both 7 ^ and Vij have the 
property 

lim 7 j (ri,r 2 )= lim 7 y(ri,r 2 ) = 0 (37) 

ri—¥oo r2—¥oo 


The color structure preserving elements of T are just 
common kinetic energy operators. 


TAi,Aj — Sij{EQ^ 


(38) 


while the elements between states with different color 
structures are given by 


1 


Tiijij — g7ii(ri324,1-4423) ~ 


2 /i 


II '' 14231 


1 

+ g(^ 0 *- ^V?324)77 (1-1324, r4423) -(39) 


B. Flux conservation 

Since T, V and g are hermitian,we have the continuity 
equation (for stationary states), 

= 0 , (40) 

or, in integral form, 

= (41) 

Oi 

where 



and the 2 factor is just our convention. 

Using the structure of T and the asymptotic behavior 
of 77 given by Eqs. (38), (39) and (37) we can easily 
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prove that non-diagonal terms on Eq. (40) do not con¬ 
tribute to the fluxes $ 0 ,. In his way, we get, 

$„ = f dSh^- • (43) 

Expanding each component as, 

( 7*') 

V’“(r) = - -Yim{0a,'^a) , (44) 


The Mqj can be chosen to form an orthonormal basis, 

(uoikoi) = > (54) 

which imposes that the parameters Aia and (pia must 
comply with the relation, 

^ ( AiaAjo^ COsijPia Pja) — d-ij . ( 55 ) 


and taking the r —>■ oo limit we obtain, 

7 ? cJn^ 

$“ = lim —. (45) 

r->oo fi^ dr 

Asymptotically u“(r) behaves as, 

u“(r) ^ sin(A:„r - ^ + (p„) + . 

(46) 

The term appears because T mixes the different chan¬ 
nels. Replacing Eq. 46 on Eq. 45, one gets. 


^ ^ . (47) 

o V a 

This is the optical theorem for our system. Considering 
for instance the case of one channel it reads. 


Replacing Eqs. (50) and (51) in Eq. (45), we are led 
to a generalization of Eq. (47), 


E 


I V' A . -f . _ A . p -f. ^ 

^ol J Jot J ict _ X '' ^oc p* p 

2z = 4" ■ 

(56) 

We note that the left side of this equation has the form, 

(57) 


_.T-TK 


2 i '■ 2i ’ 

which suggests the definition of the T matrix. 




JjOc • 


(58) 


To verify Eq. (58), we use the relation. 


Imt = |t|^ 


with 


^ = Sap , (59) 


t = 



(49) 


which can be proved from the completeness relation. 


X] IV'n)(V'n| = i , ( 60 ) 


C. T matrix 

To generalize this result, we first note that the wave- 
function can be expanded as, 

w“(r') = ^c,<(r) , (50) 

I 

with, 

M“(r) -> Aia,f^ Sili{kar - 4“ + Pia) . 

y ^cx " 

(51) 

Each uf(r) can be decomposed in two parts, 

<(r-)=<(r)+u“(r) , (52) 

where the function UQ^(r) comes from the eigenfunction 
tpg of T, and has the asymptotic limit, 

Ww(^) Aa\f^ sin{kar - + Pia) , (53) 

y '^Oc " 


of the eigenvectors of the kinetic energy T operator. With 
Eq. (59(, we calculate T^T, and indeed we prove that it 
is equal to the right side of Eq. (56) , 

TtT = ^ ’^fUja ■ (61) 

a 

Therefore, by the definition of T and the previous rela¬ 
tion, we can write relation 56 as 

ImT = TlT (62) 

which proves the unitarity of the S matrix defined by 
S = 1 -k 2tT. 


D. Identical quarks and identical antiquarks 

So far, our framework is general for four-quark systems. 
We now specialize our study to a system of two identical 
quarks and two identical antiquarks. 
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The total wavefunction must to be antissymmetric un¬ 
der the quark-exchange Pi 2 and antiquark exchange P 34 , 

Pi 2 ^' = , 

P 34 T = -T . (63) 

Given a generic wavefunction T we construct a com¬ 
pletely antissymmetric on by applying a projection 
operator 


'I'a — A/'(1 — P 12 — P 34 + Pl2-f34) '1' (64) 

In this work, our Hamiltonian is spin-independent, and 
so we ignore spin effects in the dynamics of the system. 
With this approximation, the spin part of the wavefunc¬ 
tion factorizes, and we can neglect it’s existence for ev¬ 
erything except for the symmetry of the wavefunction. 
Now, consider the function, 

4' = </>a(Pl3)«^'/3(P24)V’(l'l324)C7S (65) 

where E is the spin part of the wavefunction and ip has 
parity 


^(-r) = (-l)^’'^(i’) ( 66 ) 


We can make this function anti-symmetric by using 
Eq. (64). In this case, we obtain the function, 

Ta = Af (paiPl3)(p0{P24)tp{ri324)Cl 

_(_l)l+S'i2-7i.^^(p23)</,^(pi4).i/,(ri423)C77 

-{-l)^+^^‘^(Pa{pi4)(Pl3{p23)lp{ru23)ClI 

-^(_l)Si2-7S34-H..0^(p24)(/,^(pi3)l/,(ri324)C/ 

Ie. (67) 


We assume that E is an eigenfunction of both S '12 and 
S' 34 . This is consistent with our approximation of ne¬ 
glecting all spin related interactions, which implies that 
all spin operators commute trivially with the Hamilto¬ 
nian. Defining, 

S = [-l)Sl2 + S34 + Lr ^ 

e = (- 1 )^^^ , ( 68 ) 

where S '12 and S '34 are respectively the spins of the 12 
and 34 diquarks, we are led to a simpler expression for 

4'a, 

'I'a = A/'[4)(pi 3, p24)'*/’(l’l324)C/ + 

C4’(Pi 4, P23)V’(l’l423)C//] 51 , (69) 

with function 4) defined as, 

4>(x, y) = (pa{x)(pi3{y) + S ^a(y)(('/3(x) , (70) 

and having the property, 

(71) 


S 12 

<534 


a(-l)^’- 

J 

0 

0 

-11 

-11 

L 

0 

1 

-1 

-1 

|I/ — 1| to L + 1 

1 

0 

-11 

-1 

|I/ — 1| to L + 1 

1 

1 

-1 

-11 

L - 2 to L + 2 


Table I: Quantum numbers of the qqQQ system 


m 

Li 

n2 

L 2 

Lr 

S 

0 

0 

0 

0 

0 

-11 

0 

0 

0 

1 

1 

-1 

0 

0 

1 

0 

0 

-11 

0 

1 

0 

1 

0 

-11 


Table H: Scattering channels used for the qqQQ system. 


Since both S 12 and S 34 are conserved in the non¬ 
dynamic spin approximation, we can diagonalize our 
Hamiltonian by blocks with fixed values of these two op¬ 
erators. This is done by replacing Eq. (67) on Eq. (28). 
For this system, we obtain the scattering equation, 

fai3lp^{r)+V^fj1p^{r)+^ J d^Y'Vafi{v,Y')lp<^{Y') = 

E{ip°‘+^ J d^r'7a/3(r,r')'!/'^(r')) • (72) 

Moreover since the total orbital angular momentum 
L and parity are also conserved we can determine the 
possible "real" quantum numbers of our system. 

For instance for S '12 = S 34 = 0, we only have J = L 
states. In this case, ^ = 1 and s = (—1)'^’'. 

For S 12 = S 34 = 1, we have ^ = — 1 and still s = 
(_l)Lr total angular momentum J can range from 
|L — 2| to T -I- 2, as shown in Table I. 


E. Numerical technique 

We consider the intra-mesonic potential to be of the 
funnel type, with no other corrections. Moeover, we use 
a relativistic kinetic energy for the quarks in a meson, 

Eq = + + . (73) 

As for the scattering problem, the kinetic energy of the 
mesons is considered to be non-relativistic, but with the 
reduced mass obtained from the meson masses and not 
from the quark and antiquark masses, 

“ Mf -t Ml* ’ 

with M“ being the mass of the mesons. 

Moreover, we neglect all other relativistic effects on 
our model, such as quark-antiquark annihilation, spin- 
spin and spin-orbit interactions. Besides this, we assume 


4>(y,x) = s4-(x,y) . 
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that we are well below the threshold for decay into a 
baryon-antibaryon system. In this way we can neglected 
all the baryon-antibaryon channels, and consider just the 
meson-meson ones. 


1. Computation of meson-meson interaction 

To calculate the wave-functions of the mesons we use 
an harmonic oscillator variational basis, where each func¬ 
tion is given by, 

. (75) 

are generalized Laguerre polynomials and /? is 
a variational parameter whose value we choose as /3 = 
\ 1/3 

4m<7 \ 

3v^/ 

The Meson Hamiltonian matrix elements are calcu¬ 
lated in this basis and the matrix is diagonalized. With 
the resulting eigenfunctions, we calculate the local 
and non-local Vij parts of the meson-meson potential, 
defined by Eq. (36) and the non-local parts of the met¬ 
ric matrix ^ij defined by Eq. (31). These 7 (non¬ 
local) or 8 (local) dimensional integrals are calculated 
using Lebedev Quadrature in the angular coordinates 
and the Gauss-Laguerre Quadrature for the radial coor¬ 
dinates. We compute these multidimensional integrals 
using GPUs and the GUDA language. The resulting 
functions depend on one radial coordinate in the local 
case and on two radial coordinates on the non-local case. 
So, in both cases, the functions are evaluated on a 9- 
dimensional space. 

2. Meson-meson scattering 

We must solve Eq. 28, to obtain the asymptotic behav¬ 
ior of the wavefunction and compute the T matrix. This 
is done by discretizing the radial scattering coordinates 
on a finite box. 

We first solve the "free" equation 

f ^-0 = Eg-^o , (76) 

This is not as simple to solve as it may seem, since we 
cannot solve this equation separately for each channel, 
given the form of T and g matrices, which link all the 
channels. This is more cumbersome than what would 
happen in systems were there is only one type of quark 
combinations in mesons [2], i. e. always {giq‘i){q 2 qi) and 
never {qiqA){q 2 qz)- 

For each energy, we generate different wavefunctions 
by varying the boundary conditions of the open chan¬ 
nels. The generated solutions are not orthogonal, how¬ 
ever, and must be orthogonalized. Special care must be 


taken when doing this, since, what we want is the asymp¬ 
totic behavior of the functions. We can’t consider the in¬ 
ternal product of two wavefunctions to be just given by 
the values on the finite box used for numerics, 

{u\v) = Y^u*v, , (77) 

I 

because, although the two functions u and v are orthogo¬ 
nal in the box, the continuations of them beyond the box 
it and v aren’t necessarily orthogonal, 

{u\v) ^ 0 . (78) 

Instead, we must utilize an inner product that takes 
into account the behavior of the wavefunction continued 
outside the box. To achieve this, we consider the inner 
product of two functions to be given by 

^ ^ C0s((/?2ct i^ja) (^^) 

Oi 

in accordance with Eq. (55). The parameters A^q and 
ipia are given by Eq. (46) and are computed from the 
value of the wavefunction components at the boundary 
of the box. 

Having generated our basis of Nopen (the number of 
open channels) eigenfunctions of T for a given energy, 
we orthogonalize the basis with the Gram-Schmidt pro¬ 
cedure, applied with the inner product of Eq. 79. 

Then, for each function of the orthogonalized basis, we 
solve the scattering Eq. (28) considering 4'^ = 'I'oi + Xi- 
In this way the equation becomes 

(f + V)x^ = Egxi - V^o ^. ( 80 ) 

We solve this equation by discretizing it, using suitable 
boundary conditions for Xi- The values of fia are then 
calculated, from the behavior of Xia at the boundary of 
the box. The T matrix is evaluated utilizing Eq. (58). 

3. Finding bound states 

To find the bound states, of our system we don’t need 
to calculate the T matrix. We can directly solve Eq. (28) 
and search for states with energy smaller than the first 
threshold. 

The simplest method consists in discretizing Eq. (28) 
and solving it with Dirichlet boundary conditions. How¬ 
ever, this procedure only works as long as the spatial 
extent of the bound state is much smaller than the box 
on which we are solving the equation. When this doesn’t 
happen, the wavefunction can become highly distorted 
by the forced boundary conditions and the energy can 
even be moved above the threshold, hiding the nature 
of the state. In this work, as will be seen next, we can 
have bound states with a very small biding energy. For 
this states we would have to use a very large box to ob¬ 
tain the energy and wavefunctions accurately. However, 
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this is not convenient. Thus, instead of using Dirichlet 
boundary conditions, we consider boundary conditions 
that depend on the energy and match the expected be¬ 
havior of the wavefunctions at large distances. 

With this method, we have to solve the matrix equa¬ 
tion. 


[H + F{E)]u = Egu , (81) 

where F{E) is a matrix that fixes the boundary condi¬ 
tions. This is not a simple eigenvalue problem as the 
matrix depends on the eigenvalue itself. 

To solve this equation, we consider it as a root finding 
problem det[iJ -f F{E) — Eg\ = 0 and use the Newton’s 
method. 


_ ^(n) _ 




(82) 


to solve it. 

Applying this method to Eq. (81), we obtain the iter¬ 
ation procedure 



Figure 9: Wavefunctions of the bound states we find in 
the xxbb system. 


Tpin+l) _ pin) _ 

Fr[{H + B{E)-Eg)-^{B'{E)-g)] ’ 

(83) 

and we use it to calculate the bound-state energy. The 
wavefunction can then be calculated by solving Eq. (81) 
as a simple linear system. 

With this method, we find the bound states for the 
xxbb detailed in Section V. 


mq(GeV) 

R(MeV) 

R = ^(fm) 

1.30 

-0.02 

18.8 

1.00 

-0.95 

2.60 

0.70 

-7.91 

0.92 

0.40 

-48.54 

0.38 


Table III: Bound states in the system qqbb with L = 0 
and P = +. We use nib = i.lGeV. 


4- Finding Resonances 

To find our resonance, we extend the T matrix to the 
complex energy plane. Resonances correspond to com¬ 
plex poles of the S (or T) matrix. So, to find resonances, 
we search for zeroes of the quantity 

y[E) = l/Tr[T] . (84) 

Note that, on a pole, the T matrix is divergent and so 
is it’s trace, therefore y{E) = 0. Thus, to find a pole we 
just apply the Newton’s method to Eq. (84). 

V. RESULTS 
A. Bound states 

We apply the described method to find bound states 
in the the qqbb system. We study different values of the 
light quark mass, for angular momentum L — 0 and we 
are able to find bound states for quark masses ranging 
from niq = 0.4 to niq = 1.3 GeV. We are able to find 
bound states with s(—1)^'’ = -1-1, ^ = -1-1, listed in Table 
III. The boundstate is barely bound for = 1.3 GeV, 
the Gharm mass, and becomes more strongly bound as 


niq decreases. Wavefunctions for the first (see Table II) 
dominant channel are given on Fig. 9. 

Now, comparing Table I, these boundstates correspond 
to having 512 = ^34 = 0, and hence to L = 0 and J = 0. 
If the lightest quarks are u and d quarks, we also have 
to consider the isospin symmetry, contributing an addi¬ 
tional (— 1 )^ 12 +^ factor for the P 12 symmetry of the wave- 
function. In this way, the previous results are unchanged 
if /12 = 1. However if / = 0, the flavor wavefunction 
is anti-symmetric and so the spin wavefunction has to 
change it’s symmetry in order for the total wavefunction 
to be completely anti-symmetric. So, for I 12 = 0, we 
have also 5 i 2 = 1, and hence 5=1 and J = 1. 

To summarize, we obtain tetraquark bound states on 
the qqbb system, with quantum numbers O’*" for s, c or 
b quarks, or light quarks with 1 12 = 1. For light quarks 
with /12 = 0 , we obtain bound states with quantum num¬ 
bers 1 +. 

We also tried to find bound states for the qqcc system, 
but we were unable to find them when the lightest quarks 
have constituent masses equal or larger than the ones of 
light quarks niq > 400MeV. 
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^open 

A(GeV) 

En 

En + 1 

1 

12.6195 - 0.1248i 
12.6797- 0.16061 

12.5856 

12.974 

2 

12.9824 - 0.0078i 
12.9983 - 0.0179i 

12.974 

13.110 

3 

13.1385 - 0.07011 
13.2121 - 0.08351 

13.110 

13.3624 


Table IV: Resonances in the O’*' ccbb system, with 
^ = +1, for different numbers of open channels Nopen- 
En and En+i are the energies of the two thresholds. 


B. Resonances 

With the effective meson-meson potential, we more¬ 
over compute the T Matrix and search for poles in the S 
matrix, for the qqbb system. 

Searching for poles of the S matrix in the complex en¬ 
ergy plane for the ccbb system with L = 0, P = + and 
^ = -1-1, we find several resonances between all the four 
thresholds considered, as listed in Table IV and showed 
in Fig. 11. The pole positions correspond to the two nar¬ 
rowest resonances of each threshold interval. The narrow¬ 
est of all resonances appear between the opening of sec¬ 
ond and third channels (that Ngpen = 2 open channels). 
Of the two narrowest resonances we find in this interval, 
one has a width of 16MeV and the other 36MeV. For 
the resonances found on the other two intervals, the nar¬ 
rowest of them all has a width of 140 MeV, much larger 
than the ones found for Nopen = 2. We remember that 
the second channel is the only of the four considered ones 
that has an orbital angular momentum between the final 
state mesons. It is possible that the centrifugal barrier 
between the mesons is what is slowing the decay of these 
tetraquark resonances, as discussed in Ref. [2]. 

In Table V, we list the two narrowest resonances be¬ 
tween the second and third threshold, for different iriq 
in the O’*' ccbb system. For all values of Trig, the nar¬ 
rowest of the two resonances has also the smallest Re E. 
It is interesting to note that, the width of these two 
resonances are inversely related. When one becomes 
larger the other decreases, with varying m^. Depend¬ 
ing on the mass of the lightest quark, the width of these 
resonances can be smaller than 1 MeV or larger than 
200 MeV. Both extreme cases happen for an unphysical 
mass niq = 700 MeV. 

Broader resonances are also be found. In Fig. 11 we 
show a plot of the phase of Tr T for complex energies 
with real parts between the second and third thresholds 
and negative imaginary parts. Four resonances are there 
observed, all of them with widths smaller than 120 MeV. 
It is interesting to note the position of the four resonances 
in the complex plane forming a straight line. 

Doing the same study for ^ = —1, we also find sev¬ 
eral resonances. In Table VI, we compare these reso¬ 
nances for the ones with ^ = -|-1 for Nopen = 2 and 
TUq = me- There, we can see that, in contrast to what 


mq{GeV) 

A(GeV) 

E 2 

E 3 

1.30 

12.9824 - 0.00781 
12.9983 - 0.01791 

12.974 

13.110 

1.00 

12.4918 - 0.00481 
12.5047- 0.01921 

12.4849 

12.6241 

0.70 

12.0329 - 2 X 10"®i 
12.0496 - 0.02151 

12.0348 

12.1762 

0.40 

11.6163 - 0.00071 
11.6661 - 0.01711 

11.655 

11.796 


Table V: Resonances in the 0+ qqbb system, with 
^ = -|-1, with varying quark mass m^. E 2 and E 3 are 
the energy of the thresholds of the second and third 
channels. 



12.98 13 13.02 13.04 13.06 13.08 13.1 


Figure 10: Plot of Tr T, between the second and third 
thresholds {Nopen = 2), for the cd)b system with L = 0, 
P = + and ^ = -1-1. Note the existence of four 
resonances with r/2 < 60 MeV. 


happens for the bound states, not only resonances exist 
for 5 = — 1, but also have a behavior similar to that of 
^ = -fl. The real parts of the pole positions are very 
close and the imaginary parts are of the same order of 
magnitude. This result can be understood if we note 
that, what makes the results differ for the two values of 
^ is the presence of the non-local potential. When we 
are above the threshold, the wave-function becomes os¬ 
cillatory and the convolution of the wave-function with 
the non-local potential, vanishes if the potential varies 
slowly. In this way, only the local potential becomes im¬ 
portant well above the threshold and, consequently, the 
behavior of two resonances becomes similar. 

For the qqcc we are also able to find some resonances. 
The narrowest of them were between the opening of the 
second and third channels, similarly to what happens on 
the system. A study of these resonances for Nopen = 
2 and varying iriq is presented on Table VII. We find that 
the most stable of the two resonances has a width that 
does not vary much with niq , having always a value in the 
range of 35 — 41 MeV. The width of the second in contrast 
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^open 

? = +l 

A(GeV) 

? = -l 

1 

12.6195 - 0.1248i 

12.6069 - 0.0968i 


12.6797- 0.1606i 

12.6659 - 0.1287i 

2 

12.9824 - 0.0078i 

12.9823 - 0.0106i 


12.9983 - 0.0179i 

13.0070 - 0.0256i 

3 

13.1385 - 0.0701i 

13.1272 - 0.0341i 


13.2121 - 0.0835i 

13.1532 - 0.0682i 


Table VI: Comparison of pole positions found for 
^ = +1 and ^ = — 1 for the 0 + qqbb system. 


mq{GeY) 

A 

B 

C 

1.60 

-12.4 

- 

- 

1.30 

-13.3 

- 

-0.02 

1.00 

-16.0 

-0.2 

-0.95 

0.70 

-25.1 

-5.7 

-7.91 

0.40 

-61.2 

-41.3 

-48.54 


Table VIII: Binding energies (in MeV) for the 0+ qqbb 
system, ^ = + 1 , with three different models: A - 
Color-independent simple flip-flop, B - Color-dependent 
simple flip-flop, C - Color-dependent triple flip-flop. 


mq{GeV) 

A(GeV) 

E2 

E3 

1.30 

6.5036 - 0.0182i 
6.5427 - 0.0380i 

6.48644 

6.65375 

1.00 

5.9895 - 0.0194i 
6.0310 - 0.0413i 

5.97087 

6.14272 

0.70 

5.5113 - 0.0203i 
5.5547 - 0.0459i 

5.49134 

5.66693 

0.40 

5.0990 - 0.0184i 
5.1433 - 0.0541i 

5.07963 

5.25557 


Table VII: Resonances for the 0^ qqcc system, with 
^ = -|-1 and varying quark mass m^. 



Figure 11: Plot of Tr T, between the second and third 
thresholds {Nopen = 2 ), for the cccc system with L = 0 , 
P = + and ^ = -1-1. Note the existence of three 
resonances with r/2 < 70 MeV. 


increases monotonically from 76 MeV to 108 MeV as the 
mass TUq decreases from 1.3 GeV to 0.4 GeV. 


C. Comparison with other models 

We compare our model (C) with other two models for 
this system, defined in Section HID, namely one sim¬ 
ilar to our triple string flip-flop potential, but without 
a tetraquark sector in the potential( B) and the simple 
colorless flip-flop model (A). The results of these three 


models are compared in Table VIII. As one sees there, 
when we remove the tetraquark sector from our model, 
we still have bound states, for sufficiently small quark 
masses. For the color independent simple flip-flop model, 
the binding is even bigger than in our model, even though 
the tetraquark sector is not present at all. 

Therefore, we see that the presence of a tetraquark sec¬ 
tor in the potential is not necessary for the existence of 
bound states in the exotic tetraquark sector. This result 
has been observed by other authors, namely [65, 66 , 70] 
where bound states are found for a string flip-flop poten¬ 
tial without any tetraquark configuration considered. 

All these bound state results are for s(—1)^'’ = -|-1 and 
^ = -1-1 (see Table I), which correspond to S '12 = S '34 = 0, 
and so = O’*' for s, c or b quarks or, if 1 12 = 1, 5 'i 2 = 1 
and 5'34 = 0. However, for /12 = 0, exchange symmetry 
imposes S'! 2 = 1 , which gives quantum numbers 1 +. 


VI. DISCUSSION AND CONCLUSION 

In this work we briefly review recent experimental and 
lattice QCD results on tetraquarks. We extend the ex¬ 
isting techniques to solve tetraquarks in quark models, 
fully unitarizing for the first time the triple string flip- 
flop potential to study boundstates and resonances in 
light-light-antiheavy-antiheavy systems qqQQ. Consis¬ 
tently with the previous computations with simpler flip- 
flop potentials, we And several tetraquark boundstates 
and resonances [65, 66 , 70[. 

Comparing with recent lattice QCD results, [15, 71, 82, 
83[, which And boundstates but are not yet able to ad¬ 
dress resonances, we And a qualitative dynamical agree¬ 
ment in the sense that biding is favored when the light 
quark q gets lighter and the heavy antiquark Q gets heav¬ 
ier. 

However, concerning the symmetry and quantum num¬ 
bers, our results contradict the very recent lattice QCD 
results in Refs. [15, 71, 82, 83[. In lattice QCD, attrac¬ 
tion is only found in scalar isosinglet and vector isotriplet 
channels, while here we only And bound states (and hence 
the maximum attraction) 512 = 0 and J = 1 (scalar 
isotriplet) and 5 i 2 = 1 and thus to J = 0 (vector isosin¬ 
glet) channels. Also note that the lattice results are con¬ 
sistent with the theoretical predictions of bound qqQQ 

























16 


systems with Q heavy enough [3-13], who imply the 
heavy anti-diquark QQ color wavefunction is a triplet 
3. 

Notice we verified that the presence of the tetraquark 
sector in the flip-flop potential is not required for our 
tetraquark bound states. Indeed, we note that the so¬ 
lutions with ^ = -|-1 are mostly color symmetric, while 
those with ^ = — 1 is mostly color anti-symmetric. This 
means our bound states are mostly color symmetric, con- 
trarily to what one would expect, but consistent with 
the fact that removing the tetraquark sector from the 
potential does not destroy all the bound states. As for 
the lattice results, they are indeed color anti-symmetric 
as one would expect. The attractive channels have the 
same symmetry for spin and isospin as they are either 
scalar isosinglet (S '!2 = 0 and / = 0) or vector isotriplet 
(S '12 = 1 and J = 1) and so, as the space symmetry is 
even, the color symmetry has to be negative for the total 
wavefunction to be anti-symmetric. 

Let us analyze the mechanism why we obtain attrac¬ 
tion for the color symmetric case and repulsion in the 
color anti-symmetric case. For the qqQQ system, if we 
consider only two channels (related by quark and anti¬ 
quark exchange), the scattering equation becomes 


T 2 p{r) + V^^{r) + ^ J (fr'v{r,r')'>p{r') = 
E{ip{r)+^ J ,) (85) 


calculated from Eq. 36. We see that it depends on E/j/ — 
. When the ground state is vq = Vj, we have 
Viji = ^Vi and so, this becomes . Since Vi is the 

ground state, we have have Vj < Vn and, therefore 


Vi - Vii 

6 


< 0 


( 86 ) 


Since the function $ in the integrand is node-less (be¬ 
cause it is the ground state). 


‘J’(Pl3) P24)*‘h(pi4, P 23 ) > 0 . (87) 


Therefore, in this limit we have 


u(r, r') < 0 


( 88 ) 


for r' —>■ 00 and fixed r. The same result occurs if we 
exchange r and r'. Then, in Eq. (85) the non-local po¬ 
tential will be mostly attractive for ^ = -1-1 and mostly 
repulsive for ^ = — 1. This confirms our results. 

To conclude we develop fully unitarized techniques to 
study tetraquarks with quark models, and to search for 
boundstates and resonances. Asymptotically, the four 
quark system with a string flip-flop potential reduces to 
coupled two-body meson-meson systems, with non-local 
potentials that vanish at long distances, thus solving the 
Van der Waals problem. However, we And that the string 
flip-flop potential still remains attractive enough to pro¬ 
duce qqQQ boundstates with quantum numbers not ob¬ 
served in recent lattice QCD computations. 

It should be noted that different solutions to remove 
the excessive attraction exist [84]. For instance, attrac¬ 
tion may change if the non-local potentials change sign 
in the region where ri ~ r 2 , and the presence of the gafi 
operator cancels in part this effect. We expect that our 
results will motivate future studies of different solutions 
to this excessive attraction. 

Marco Cardoso is supported by FCT under the con¬ 
tract SFRH/BPD/73140/2010. 
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